###################################################################################
# This source code file reproduces the analysis of
# "The Most Liberal Senator" in No News is News
# It calls four different files:
# 	national_journal_2003_votes_bush.txt  # Raw data from the National Journal
#	OnlyRelevantSenatorsIRT.RData		  # Posterior distribution IRT model
#	OnlyRelevantSenatorsNEG.RData		  # Posterior distribution, abstentions as NAYs
#	OnlyRelevantSenators.RData			  # Competing principals, informative priors
#	OnlyRelevantSenatorsFLATPRIORS.RData  # Competing principals, FLAT priors
###################################################################################


# Require the following packages
library (runjags)
library (mcmcplots)
library (random)
library (multicore)

# Load votes
njkey <- read.table("national_journal_2003_votes_bush.txt", header=TRUE, sep="\t", na.strings=9)
poo <- apply(is.na(njkey[,-c(1:3)]),1,sum)/dim(njkey[,-c(1:3)])[2]       #   Calculate Percent Missing

names <- strsplit (as.character(njkey$name), split=c(" ","  "))
new.name <- c()
for (i in 1:101) {
  x <- names[[i]][1]
  new.name <- c(new.name, x)
}

party.ab <- ifelse (njkey$party==200, "R", ifelse (njkey$party==300, "D", "D"))
senate.id <- paste (new.name, " (", party.ab, "-", njkey$statenm, ")", sep="")
Roll.Call <- njkey[-1,-c(1:3)]
for (i in 1:nrow(Roll.Call)){
  for (j in 1:ncol(Roll.Call)){
    Roll.Call[i,j] <- ifelse (Roll.Call[i,j]==6 | Roll.Call[i,j]==8, 0, Roll.Call[i,j])}}

Roll.Call <- as.matrix (Roll.Call)
dim (Roll.Call[party.ab[-1]=="R", ])
dim (Roll.Call[party.ab[-1]=="D", ])
Rep.rolls <- colSums( Roll.Call[party.ab[-1]=="R", ], na.rm=T)/51
Dem.rolls <- colSums( Roll.Call[party.ab[-1]=="D", ], na.rm=T)/49

cbind (Rep.rolls, Dem.rolls)
fix.bill <- c(which (colnames(Roll.Call)=="X115"), which (colnames(Roll.Call)=="X180"))
fix.bill <- sort (fix.bill)

# Identify abstainers within Democratic party
kerry   <- which(senate.id[-1]=="KERRY (D-MA)")
graham  <- which(senate.id[-1]=="GRAHAM (D-FL)")
edwards <- which(senate.id[-1]=="EDWARDS (D-NC)")
lieber  <- which(senate.id[-1]=="LIEBERMAN (D-CT)")

senatorsRunning <- sort ( c(kerry, graham, edwards, lieber)) 
new.order <- c(seq (1:100)[-senatorsRunning], senatorsRunning)

rownames (Roll.Call) <- new.name[-1]
# Party.membership variables are wrong.  There are 101, when we should have 98
party.membership <- round (njkey$party/100)  
party.membership <- ifelse (party.membership==3,2,party.membership)
# Adding the next two lines corrects the mistake with party.membership
party.membership.no.bush <- party.membership[-1]
party.membership.no.bush <- party.membership.no.bush[new.order]
lead.L  <- Roll.Call["REID",]     # Reid is DP whip
lead.R  <- Roll.Call["NICKLES",]  # Nickles is GOP whip
lead <- rbind(lead.L, lead.R)

# Crucial step here: reorder Roll.Call
Roll.Call.reorder <- Roll.Call[new.order,]
rc.nor <- ifelse(is.na(Roll.Call.reorder), NA
		, ifelse(Roll.Call.reorder==0, 0, 1))
rc.abs <- ifelse(is.na(Roll.Call.reorder), 2
		, ifelse(Roll.Call.reorder==0,1,3))  # code NA as 2, Nay as 1, Aye as 3

rc <- rbind (rc.nor[-c(97:100),], rc.abs[c(97:100),])
rc2 <- rc[-c(which (rownames(rc)=="REID"), which (rownames(rc)=="NICKLES")),]
new.senate.id <- senate.id[-1][new.order][-c(which (rownames(rc)=="REID"), which (rownames(rc)=="NICKLES"))]

party.membership.no.whips <- party.membership.no.bush[-c(which (rownames(rc)=="REID"), which (rownames(rc)=="NICKLES"))]
n.legislators <- nrow (rc2)
n.item <- ncol (rc2)


########################################################################################
# This is the regular IRT model, run for comparison, using exact same priors
########################################################################################
irt="model { 
  for (i in 1:n.legislators) {
  for (j in 1:n.item) {
    rc[i,j] ~ dbern(mu[i,j]);
    probit(mu[i,j]) <- beta[j]*theta[i] - alpha[j];
    }
  }
# PRIORS
for(j in 1:n.item) { alpha[j] ~ dnorm(0, 0.25); }
for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0, 0.25); }
    beta[fix.bill[1]] <- 2
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
    beta[fix.bill[2]] <- -2
for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
for(i in 1:n.legislators)  { theta[i] ~ dnorm(0, 1); }
}"

y <- rc.nor[-c(which (rownames(rc)=="REID"), which (rownames(rc)=="NICKLES")),]
n.legislators <- nrow (y)
n.item <- ncol (y)

rc.against <- matrix (NA, nrow=nrow(y), ncol=ncol(y))
for (i in 1:nrow(y)){
  for (j in 1:ncol(y)){
    rc.against[i,j]   <- ifelse (is.na(y[i,j]), abs(lead[party.membership.no.whips[i],j]-1), y[i,j])  
  }
}
inits.theta <- princomp (rc.against)$scores[,1]

irt.parameters = c("theta", "alpha", "beta", "deviance")
irt.inits <- function() {
  dump.format(
    list(
      theta = inits.theta + rnorm(n.legislators)
      , alpha = rnorm(n.item)
      , beta = c(rnorm(fix.bill[1]-1), NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}
irt.data <- dump.format(list(rc=y, n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill))

# Uncomment the following lines to run the model, with three chains distributed among three cores
# This process takes a very long time
# The posterior distribution of parameters appears in file:
# "OnlyRelevantSenatorsIRT.RData"
# system.time(
#   out3 <- mclapply(1:3, function(x) {
# 		jags.irt <- try(run.jags(   
# 			model=irt,
# 			monitor=irt.parameters,
# 			n.chains=1,
# 			data=irt.data,
# 			inits=irt.inits(),
# 			thin=50, burnin=25000, sample=50000,
# # 			      thin=5, burnin=200, sample=200,
# 			check.conv=FALSE, plots=FALSE
# 		))
# 		if(inherits(jags.irt,"try-error")) {return()}
# 		return(jags.irt)
# 	}, mc.cores=3 )
# )
# 
# save (out3, file="OnlyRelevantSenatorsIRT.RData")  # informative priors


########################################################################################
# NA's as votes against the whip.
########################################################################################
n.legislators <- nrow (rc.against)
n.item <- ncol (rc.against)

neg.parameters = c("theta", "alpha", "beta", "deviance")
neg.inits <- function() {
  dump.format(
    list(
    	theta = inits.theta + rnorm(n.legislators)
	  , alpha = rnorm(n.item)
      , beta = c(rnorm(fix.bill[1]-1), NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
  )
}
neg.data.against <- dump.format(list(rc=rc.against, n.legislators=n.legislators, n.item=n.item, fix.bill=fix.bill))

# Uncomment the following lines to run the model, with three chains distributed among three cores
# This process takes a very long time
# The posterior distribution of parameters appears in file:
# "OnlyRelevantSenatorsNEG.RData"
# system.time(
#   out2 <- mclapply(1:3, function(x) {
#     jags.neg <- try(run.jags(   
#       model=irt,
#       monitor=neg.parameters,
#       n.chains=1,
#       data=neg.data.against,
#       inits=neg.inits(),
#       thin=50, burnin=25000, sample=50000,
# #       			      thin=5, burnin=200, sample=200,
#       check.conv=FALSE, plots=FALSE
#     ))
#     if(inherits(jags.neg,"try-error")) {return()}
#     return(jags.neg)
#   }, mc.cores=3 )
# )
# 
# save (out2, file="OnlyRelevantSenatorsNEG.RData")  # informative priors





########################################################################################
# This is the new model with abstention process only for four democrats
########################################################################################
compete.bis="model { 
  for (i in 1:(n.legislators-4)) {
  	for (j in 1:n.item) {
    	y[i,j] ~ dbern(mu[i,j]);
    	probit(mu[i,j]) <- beta[j]*theta[i] - alpha[j];
    }
  }
  for (i in (n.legislators-3):n.legislators) {
    for (j in 1:n.item) {
      probit(CJR[i,j]) <- beta[j]*theta[i] - alpha[j];
      y[i,j] ~ dcat(Q[i,j,1:3]);         # code Nay as 1, NA as 2, Aye as 3
      Q[i,j,1] <- pow(probVoteDisagree*(1-CJR[i,j]), p[i,j]) * pow(probVoteAgree*(1-CJR[i,j]), (1-p[i,j]));
      Q[i,j,2] <- 1 - Q[i,j,1] - Q[i,j,3];
      Q[i,j,3] <- pow(probVoteAgree*CJR[i,j], p[i,j]) * pow(probVoteDisagree*CJR[i,j], (1 - p[i,j]));
      p[i,j] <- lead[party.membership[i],j];
    } 
  }
# priors
#   probVoteDisagree ~ dbeta(6,28)   # This is a specialized prior
#   probVoteAgree ~ dbeta(16,35)     # This is a specialized prior
  probVoteDisagree ~ dbeta(1,1)  # This is a flat prior
  probVoteAgree ~ dbeta(1,1)     # This is a flat prior
for(j in 1:n.item) { alpha[j] ~ dnorm(0, 0.25); }
for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0, 0.25); }
    beta[fix.bill[1]] <- 2
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0, 0.25); }
    beta[fix.bill[2]] <- -2
for(j in (fix.bill[2]+1):n.item) { beta[j]  ~ dnorm(0, 0.25); }
for(i in 1:n.legislators)  { theta[i] ~ dnorm(0, 1); }
}"

compete.parameters = c("theta", "alpha", "beta", "probVoteAgree", "probVoteDisagree", "deviance")

compete.data <- dump.format(list(y=rc2
								 , lead=lead
								 , party.membership=party.membership.no.whips
								 , n.legislators=n.legislators
								 , n.item=n.item
								 , fix.bill=fix.bill))

compete.inits <- function() {
	dump.format(
		list(
			theta = rnorm(n.legislators)
			, alpha = rnorm(n.item)
			, beta = c(rnorm(fix.bill[1]-1), NA, rnorm(fix.bill[2]-fix.bill[1]-1), NA, rnorm(n.item-fix.bill[2]))
			, probVoteDisagree = runif (1,0,0.3)
			, probVoteAgree = runif (1,0.3,0.6)
			,.RNG.name="base::Wichmann-Hill"
			,.RNG.seed=randomNumbers(n = 1, min = 1, max = 1e+04,col=1))
	)
}


# Uncomment the following lines to run the model, with three chains distributed among three cores
# This process takes a very long time
# The posterior distribution of parameters appears in file:
# "OnlyRelevantSenators.RData" and
# "OnlyRelevantSenatorsFLATPRIORS.RData"

# system.time(
# 	out4 <- mclapply(1:3, function(x) {
# 		jags.compete.new.priors <- try(run.jags(   
# 			model=compete.bis,
# 			monitor=compete.parameters,
# 			n.chains=1,
# 			data=compete.data,
# 			inits=compete.inits(),
# 			thin=50, burnin=25000, sample=50000,
# # 			      thin=5, burnin=200, sample=200,
# 			check.conv=FALSE, plots=FALSE
# 		))
# 		if(inherits(jags.compete.new.priors,"try-error")) {return()}
# 		return(jags.compete.new.priors)
# 	}, mc.cores=3 )
# )

# save (out4, file="OnlyRelevantSenators.RData")  # informative priors
# save (out4, file="OnlyRelevantSenatorsFLATPRIORS.RData")  # flat priors

#  STOP HERE, then reload objects
load ("OnlyRelevantSenatorsIRT.RData")
anchors <- c(grep(paste("beta[",fix.bill[1],"]",sep=""), colnames (out3[[1]][[1]][[1]]), fixed=TRUE)
			 , grep(paste("beta[",fix.bill[2],"]",sep=""), colnames (out3[[1]][[1]][[1]]), fixed=TRUE))

chain1 <- out3[[1]][[1]][[1]][,-anchors]  # eliminate fixed beta parameters
chain2 <- out3[[2]][[1]][[1]][,-anchors]  # and spike priors on probs of voting
chain3 <- out3[[3]][[1]][[1]][,-anchors]

# Convergence tests
chainsIRT <- mcmc.list(list (chain1, chain2, chain3))  
gelman.diag (chainsIRT)
# Traceplots of theta parameters (will return an html file)
mcmcplot (chainsIRT,    
		  , filename="traceplots2params"   # <= Change name here
		  , random=10)
rm (out3, anchors, chain1, chain2, chain3)


load ("OnlyRelevantSenatorsFLATPRIORS.RData")
anchors <- c(grep(paste("beta[",fix.bill[1],"]",sep=""), colnames (out4[[1]][[1]][[1]]), fixed=TRUE)
			 , grep(paste("beta[",fix.bill[2],"]",sep=""), colnames (out4[[1]][[1]][[1]]), fixed=TRUE))

chain1 <- out4[[1]][[1]][[1]][,-anchors]  # eliminate fixed beta parameters
chain2 <- out4[[2]][[1]][[1]][,-anchors]  # and spike priors on probs of voting
chain3 <- out4[[3]][[1]][[1]][,-anchors]

# Convergence tests
chainsFlat <- mcmc.list(list (chain1, chain2, chain3))  
gelman.diag (chainsFlat)
# Traceplots of theta parameters  (will return an html file)
mcmcplot (chainsFlat,  
		  , filename="traceplots2params"   # <= Change name here
		  , random=10)
rm (out4, anchors, chain1, chain2, chain3)


load ("OnlyRelevantSenators.RData")
anchors <- c(grep(paste("beta[",fix.bill[1],"]",sep=""), colnames (out4[[1]][[1]][[1]]), fixed=TRUE)
			 , grep(paste("beta[",fix.bill[2],"]",sep=""), colnames (out4[[1]][[1]][[1]]), fixed=TRUE))
chain1 <- out4[[1]][[1]][[1]][,-anchors]  # eliminate fixed beta parameters
chain2 <- out4[[2]][[1]][[1]][,-anchors]  # and spike priors on probs of voting
chain3 <- out4[[3]][[1]][[1]][,-anchors]
chainsConv <- mcmc.list(list (chain1, chain2, chain3))  

# Convergence tests
gelman.diag (chainsConv)
# Traceplots of theta parameters (will return an html file)
mcmcplot (chainsConv, 
          , filename="traceplots2params"   # <= Change name here
          , random=10)
rm (out4, anchors, chain1, chain2, chain3)


load ("OnlyRelevantSenatorsNEG.RData")
anchors <- c(grep(paste("beta[",fix.bill[1],"]",sep=""), colnames (out2[[1]][[1]][[1]]), fixed=TRUE)
			 , grep(paste("beta[",fix.bill[2],"]",sep=""), colnames (out2[[1]][[1]][[1]]), fixed=TRUE))
chain1 <- out2[[1]][[1]][[1]][,-anchors]  # eliminate fixed beta parameters
chain2 <- out2[[2]][[1]][[1]][,-anchors]  # and spike priors on probs of voting
chain3 <- out2[[3]][[1]][[1]][,-anchors]
chainsNeg <- mcmc.list(list (chain1, chain2, chain3))  
gelman.diag (chainsNeg)

# Traceplots of theta parameters (will return an html file)
mcmcplot (chainsNeg,    
		  , filename="traceplots2params"   # <= Change name here
		  , random=10)
rm (out2, anchors, chain1, chain2, chain3)


# EXTRACT probabilities of voting
prob.Flat <- rbind (chainsFlat[[1]][sample(1:nrow(chainsFlat[[1]]), 200), grep("probVote", colnames(chainsFlat[[1]]))]
                    , chainsFlat[[2]][sample(1:nrow(chainsFlat[[2]]), 200), grep("probVote", colnames(chainsFlat[[2]]))]
                    , chainsFlat[[3]][sample(1:nrow(chainsFlat[[3]]), 200), grep("probVote", colnames(chainsFlat[[3]]))])

prob.Conv <- rbind (chainsConv[[1]][sample(1:nrow(chainsConv[[1]]), 200), grep("probVote", colnames(chainsConv[[1]]))]
                    , chainsConv[[2]][sample(1:nrow(chainsConv[[2]]), 200), grep("probVote", colnames(chainsConv[[2]]))]
                    , chainsConv[[3]][sample(1:nrow(chainsConv[[3]]), 200), grep("probVote", colnames(chainsConv[[3]]))])

# EXTRACT deviance parameters
Deviance.Flat <- c (chainsFlat[[1]][sample(1:nrow(chainsFlat[[1]]), 200), grep("deviance", colnames(chainsFlat[[1]]))]
                    , chainsFlat[[2]][sample(1:nrow(chainsFlat[[2]]), 200), grep("deviance", colnames(chainsFlat[[2]]))]
                    , chainsFlat[[3]][sample(1:nrow(chainsFlat[[3]]), 200), grep("deviance", colnames(chainsFlat[[3]]))])

Deviance.Conv <- c (chainsConv[[1]][sample(1:nrow(chainsConv[[1]]), 200), grep("deviance", colnames(chainsConv[[1]]))]
                    , chainsConv[[2]][sample(1:nrow(chainsConv[[2]]), 200), grep("deviance", colnames(chainsConv[[2]]))]
                    , chainsConv[[3]][sample(1:nrow(chainsConv[[3]]), 200), grep("deviance", colnames(chainsConv[[3]]))])

Deviance.IRT  <- c (chainsIRT[[1]][sample(1:nrow(chainsIRT[[1]]),200), grep("deviance", colnames(chainsIRT[[1]]))]
                    , chainsIRT[[2]][sample(1:nrow(chainsIRT[[2]]),200), grep("deviance", colnames(chainsIRT[[2]]))]
                    , chainsIRT[[3]][sample(1:nrow(chainsIRT[[3]]),200), grep("deviance", colnames(chainsIRT[[3]]))])

Deviance.Neg <- c (chainsNeg[[1]][sample(1:nrow(chainsNeg[[1]]),200), grep("deviance", colnames(chainsNeg[[1]]))]
                   , chainsNeg[[2]][sample(1:nrow(chainsNeg[[2]]),200), grep("deviance", colnames(chainsNeg[[2]]))]
                   , chainsNeg[[3]][sample(1:nrow(chainsNeg[[3]]),200), grep("deviance", colnames(chainsNeg[[3]]))])


# EXTRACT Theta parameters
Theta.Flat <- rbind (chainsFlat[[1]][sample(1:nrow(chainsFlat[[1]]), 200), grep("theta", colnames(chainsFlat[[1]]))]
                     , chainsFlat[[2]][sample(1:nrow(chainsFlat[[2]]), 200), grep("theta", colnames(chainsFlat[[2]]))]
                     , chainsFlat[[3]][sample(1:nrow(chainsFlat[[3]]), 200), grep("theta", colnames(chainsFlat[[3]]))])

Theta.Conv <- rbind (chainsConv[[1]][sample(1:nrow(chainsConv[[1]]), 200), grep("theta", colnames(chainsConv[[1]]))]
                     , chainsConv[[2]][sample(1:nrow(chainsConv[[2]]), 200), grep("theta", colnames(chainsConv[[2]]))]
                     , chainsConv[[3]][sample(1:nrow(chainsConv[[3]]), 200), grep("theta", colnames(chainsConv[[3]]))])

Theta.IRT  <- rbind (chainsIRT[[1]][sample(1:nrow(chainsIRT[[1]]),200), grep("theta", colnames(chainsIRT[[1]]))]*(-1)
                     , chainsIRT[[2]][sample(1:nrow(chainsIRT[[2]]),200), grep("theta", colnames(chainsIRT[[2]]))]*(-1)
                     , chainsIRT[[3]][sample(1:nrow(chainsIRT[[3]]),200), grep("theta", colnames(chainsIRT[[3]]))]*(-1))

Theta.Neg <- rbind (chainsNeg[[1]][sample(1:nrow(chainsNeg[[1]]),200), grep("theta", colnames(chainsNeg[[1]]))]*(-1)
                    , chainsNeg[[2]][sample(1:nrow(chainsNeg[[2]]),200), grep("theta", colnames(chainsNeg[[2]]))]*(-1)
                    , chainsNeg[[3]][sample(1:nrow(chainsNeg[[3]]),200), grep("theta", colnames(chainsNeg[[3]]))]*(-1))


# Standardize ideal points so that they can be directly comparable
Standardize.Theta <- function (Theta) {
  min.t <- min (Theta)
  max.t <- max (Theta)
  dis.t <- max.t-min.t
  res <- apply (Theta, c(1,2), function (x) { (x-min.t)/(dis.t) })
  return (res)
}


Theta.sig.Conv <- matrix (NA, ncol=3, nrow=ncol(Theta.Conv))
for (i in 1:ncol(Theta.Conv)){
  Theta.sig.Conv[i,] <- quantile (Theta.Conv[,i], prob=c(0.25,0.5,0.75))
}

Theta.sig.Flat <- matrix (NA, ncol=3, nrow=ncol(Theta.Flat))
for (i in 1:ncol(Theta.Flat)){
	Theta.sig.Flat[i,] <- quantile (Theta.Flat[,i], prob=c(0.25,0.5,0.75))
}

Theta.sig.IRT <- matrix (NA, ncol=3, nrow=ncol(Theta.IRT))
for (i in 1:ncol(Theta.IRT)){
  Theta.sig.IRT[i,] <- quantile (Theta.IRT[,i], prob=c(0.25,0.5,0.75))
}

Theta.sig.Neg <- matrix (NA, ncol=3, nrow=ncol(Theta.Neg))
for (i in 1:ncol(Theta.Neg)){
  Theta.sig.Neg[i,] <- quantile (Theta.Neg[,i], prob=c(0.25,0.5,0.75))
}

Theta.Std.Flat <- Standardize.Theta (Theta.sig.Flat)
Theta.Std.Conv <- Standardize.Theta (Theta.sig.Conv)
Theta.Std.IRT  <- Standardize.Theta (Theta.sig.IRT)
Theta.Std.Neg  <- Standardize.Theta (Theta.sig.Neg)


ord.ideal <- order (Theta.Std.IRT[,2])
abs.rate <- round (rowSums(is.na(rc.nor[-c(which (rownames(rc)=="REID"), which (rownames(rc)=="NICKLES")),]))/ncol(rc.nor), 2)

# The following code produces Figure 2
# pdf ("senateComparison2013twoPar.pdf", h=13, w=9)
par (mar=c(3,1,1,1), cex.axis=1, las=1)
plot ( c(-4,3), c(1,98), type="n", ylab="", xlab="", main="", axes=F, xlim=c(0,1.25))
axis(1, at=c(0,1.25), labels=NA)
segments (x0=Theta.Std.IRT[ord.ideal,1], x1=Theta.Std.IRT[ord.ideal,3], y0=1:98, y1=1:98, col="gray")
points ( xy.coords ( Theta.Std.IRT[ord.ideal,2], (1:98)), pch=19, col="gray")
# segments (x0=Theta.Std.Neg[ord.ideal,1], x1=Theta.Std.Neg[ord.ideal,3], y0=(1:98)+0.1, y1=(1:98)+0.1, col="black")
points ( xy.coords ( Theta.Std.Neg[ord.ideal,2], (1:98)+0.2), pch=1, col="black")
segments (x0=Theta.Std.Conv[ord.ideal,1], x1=Theta.Std.Conv[ord.ideal,3], y0=(1:98)+0.2, y1=(1:98)+0.2, col="black")
points ( xy.coords ( Theta.Std.Conv[ord.ideal,2], (1:98)+0.2), pch=19, col="black")
# segments (x0=Theta.Std.Flat[ord.ideal,1], x1=Theta.Std.Flat[ord.ideal,3], y0=(1:98)-0.1, y1=(1:98)-0.1, col="black")
# points ( xy.coords ( Theta.Std.Flat[ord.ideal,2], (1:98)-0.1), pch=19, col="black")
# abline (h=which(new.senate.id[ord.ideal]=="KERRY (D-MA)"), lty=3, col="red")
legend (x=-0, y=98, legend=c("Complete-data model","Observed-data model","Abstentions as 'against party' votes"), pch=c(19,19,1), col=c("black","gray","black"), bty="n", cex=0.9)
text ( xy.coords (Theta.Std.IRT[ord.ideal,2][-1]+0.15, seq(1,98,by=1)[-1]), label=paste (new.senate.id[ord.ideal][-1], abs.rate[ord.ideal][-1], sep="-"), cex=0.6  )
text (  xy.coords ( Theta.Std.IRT[ord.ideal,2][1]+0.43, seq(1,98,by=1)[1]), label=paste (new.senate.id[ord.ideal][1], abs.rate[ord.ideal][1], sep="-"), cex=0.6, font=2)
# text (  xy.coords ( Theta.sig.31[ord.ideal,2][c(7,13,23,37)]+0.8, seq(1,98,by=1)[c(7,13,23,37)]), label=paste (senate.id[-c(1,57,73)][ord.ideal][c(7,13,23,37)], abs.rate[ord.ideal][c(7,13,23,37)], sep="-"), cex=0.6  )
# dev.off ()


#  Calculate interquartile range for Kerry's rank
kerry <- which(new.senate.id=="KERRY (D-MA)")
kerry.rank.conv <- kerry.rank.flat <- kerry.rank.irt <- numeric ()
for (i in 1:nrow(Theta.IRT)){
  kerry.rank.irt  <- c(kerry.rank.irt, rank (Theta.IRT[i,]*(-1))[kerry])
  kerry.rank.flat <- c(kerry.rank.flat, rank (Theta.Flat[i,]*(-1))[kerry])
  kerry.rank.conv <- c(kerry.rank.conv, rank (Theta.Conv[i,]*(-1))[kerry])
}
quantile (kerry.rank.irt, prob=c(0.25,0.5,0.75))
quantile (kerry.rank.flat, prob=c(0.25,0.5,0.75))
quantile (kerry.rank.conv, prob=c(0.25,0.5,0.75))





